Dark Matter Gravitational Clustering With a Long-Range Scalar Interaction 
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We explore the possibility improving the ACDM model at megaparsec scales by introducing 
a scalar interaction that increases the mutual gravitational attraction of dark matter particles. 
Using N-body simulations, we study the spatial distribution of dark matter particles and halos. We 
measure the effect of modifications in the Newton's gravity on properties of the two-point correlation 
function, the dark matter power spectrum, the cumulative halo mass function and density probability 
distribution functions. The results look promising: the scalar interactions produce desirable features 
at megaparsec scales without spoiling the ACDM successes at larger scales. 

PACS numbers: 98.80.-k, 95.35.+d, 98.65.Dx 



I. INTRODUCTION 

In this paper we study cosmological implications of a 
scalar interaction that produces a long-range fifth force 
in the dark sector, proposed by G. Farrar , S. Gubser and 
J. Peebles [E B B 0| • The physical motivation for this 
model comes from the string theory Q . 

The cosmological motivation comes from small-scale 
difhculties of the ACDM model, which successfully passes 
almost all observational tests (see e.g. Q)- Difficulties 
appear at length scales below few megaparsecs: (1) the 
ACDM voids are less empty than the real voids; (2) the 
present accretion rate of intergalactic debris onto thin 
spiral galaxies poses a problem for current galaxy forma- 
tion paradigm; (3) merger rates at low redshifts for giant 
elliptical galaxies suggest violent accretion histories for 
their haloes at low redshifts which is in contradiction 
with observations. 

The void problem, pointed out by J. Peebles 'Sj is hotly 
debated in the literature, with arguments supporting his 
orig inal claim [^, [l^l as well as arguments to the contrary 

The late merging problem appears in simulations 
which shows that accretion onto giant elliptical galaxies 
at cluster centers continues until z = [l3| while the in- 
dependence of the color-magnitude relation of SDSS[63t 
galaxies on their environment [TJ] and the remarkable 
stability of the color-magnitude relation at z = 0.7 
is consistent with the picture of giant galaxies as island 
universes, contrary to ACDM simulations. Likewise, ac- 
cording to J. Peebles, the very existence of large galaxies 
like the Milky Way, with its spiral structure intact, sug- 
gests a lack of major mergers in recent history (see [l6| 
for thin disk dominated galaxies survival issues). There 
is observational evidence that the latest " maj or invasion" 
happened to our Galaxy 10-12 Gyr ago [171. There are 
also arguments to the contrary, claiming that Milky Way- 
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like galaxies survive late merging in N-body simulations 

For an excellent discussion of the observational situ- 
ation and comparison with the ACDM model, see Refs. 
[H, US im , and the references therein. 

Theoretical suggestions of Peebles and Gubser were 
followed by the work of Nusser et al. [l^], exploring 
the cosmological implications with N-body simulations. 
Some preliminary results on similar scalar model were 
also obtained by Rodrfguez-Meza et al. [S^, [H, . The 
Long-Range Scalar Interaction (LRSI) model started a 
debate in the literature recently, mainly focused on the 
weak-equivalence principle violation and it's impact on 
dynamic of Milky- Way satellites [H [H, [l^ . 

The work, presented here should be regarded as the 
next step in this line of research. Like Nusser et at, we 
study the two-point correlation function and the statisti- 
cal properties of the mass distribution of the dark matter 
halos. To our knowledge, for the first time in the litera- 
ture, we also study the power spectra for a set of scalar 
field parameters with comoving screening length. Anal- 
ogous model was analyzed in great detail by Grawdohl 
& Frieman nearly 20 years ago [l^ll^. However, their 
model assumed a fixed physical screening length, which 
is a fundamental difference in comparison to our model. 
We have used more particles in our simulations compared 
to those used by Nusser et al., and our resolution is bet- 
ter than their resolution; we also consider a wider range 
of scalar model parameters. As a consequence, we re- 
solve the power spectrum near the screening length char- 
acteristic scale. Our results show a clear feature in the 
power spectrum near a wave number, that is the inverse 
of the screening length. The extra power seen at higher 
wavenumbers is generated by the gravitational field, en- 
hanced by the scalar interaction. Nusser et al. could not 
see a similar feature in their correlation function because 
the screening lengths they considered were too close to 
the mean interparticle separation in their simulations. 
Apart from this important difference, we confirm their 
results. The scalar field generates lower density in voids. 
It also shifts structure formation to higher redshifts. 

This paper is organized as follows: In section [TT] we in- 
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troduce the effective gravitational potential and modified 
force law used as an approximation of the scalar field. In 
Section ITTll we describe our N-body simulations. Our re- 
sults are presented in Section IIVI A brief summary and 
discussion appears in Section |Vl 



II. THEORY 

Following Refs. [1, we consider dark matter (DM) 
particles as strings. Their dynamics is defined by con- 
ventional gravity as well as an additional attractive force, 
induced by an exchange of a massless scalar. This force is 
well represented by a Yukawa-like potential with a char- 
acteristic screening length dynamically generated by the 
presence of the light particles coupled to scalar (see Ref. 
[l| for details on the dynamical screening mechanism). 
We consider one species of strings as DM particles. The 
force between two DM particles, each of mass m, arises 
from the potential 



Gm 



with 



g{x) = l + zSe-^/'-^ 



(1) 



(2) 



Here G is Newton's constant; r and x — v/a{t) are, re- 
spectively, the particle separation vector in real space and 
comoving coordinates; t is the cosmological time; a is the 
scale factor, normalized to unity at present. 



\ . 



(3) 



Here and below the subscript "0" denotes the present 
epoch. The parameter Ts is the screening length and (3 
is a measure of the relative strength of the scalar inter- 
action compared to conventional gravity. The screening 
length Ts is constant in comoving coordinates because of 
the dynamical screening mechanism, specific to the class 
of scalar fields considered here. Accordingly, modified po- 
tential gives rise to the modified force law between DM 
particles. The new appropriate form is 



DM 



TOl • 7712 



(4) 



We can adopt this modification as a distance dependent 
correction term to ordinary Newton force law: 

FdM = FMewton ' Fs{r, /S^rs), (5) 

where Fg characterize deviations from standard gravity: 



F,(r,/?,r,) = H-/3( l + -)e- 



(6) 



For /3 = or r ^ Ts we have — > 1, thus we recover 
standard Newtonian force law. 



Switching from the discrete particle picture to fluid 
dynamics, we will now introduce the dark matter density 
field, given by the expression 



p(x,i) = {p) (l + S), 



(7) 



where {p{t)) is the ensemble average of the dark matter 
density at time t, and 5(x, t) describes local deviations 
from homogeneity. The structure formation is driven 
only by the spatially fluctuating part of the gravitational 
potential, 0(x, t), induced by the density fluctuation field 
<5, 



G(p)a2 



1 — 37rf (1^-^ I) 



The Fourier transform of this equation is 



2a 



1 



where 



and 



<>k = (2^) 



(27r) 



-2/3 



-2/3 



f + {krsY 



-«k-x j3 



d^x 



d^x 



(8) 



(9) 



(10) 



(11) 



are the Fourier transforms of 4>{x) and (5(x), respectively; 
k is the comoving wavevector, and the quantities 



(12) 



and Hq are, respectively, the present values of the di- 
mensionless mean dark matter density and the Hubble 
parameter. From now on, we will also use the symbols 
h and il,A, denoting the dimensionless Hq, expressed in 
units of 100 km s~'^Mpc~^ , and the cosmological constant 
contribution to the present mean density. 

Note that when P = oi x rg, equations ^ and 
(H]) become identical with conventional gravity ^30|. In- 
deed, consider the fractional deviation from the Newto- 
nian gravitational potential. 



(13) 



Throughout this paper, the label 'Newton', and the sub- 
script 'N' refer to the ACDM cosmology with the con- 
ventional Newtonian gravity. Equation ^ gives 
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1 + {krs)- 



(14) 



This is the Fourier image of the spatial decline of the 
Yukawa potential: in the limit kr^ — > 0, the above ex- 
pression vanishes. In the opposite limit, the fractional 
deviation from Newtonian gravity reaches a finite value. 



^(l>k/<t)k P for fcr^ > 1, 
A(^k/0k for /cr^ < 1 . 



(15) 
(16) 
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We consider values of (3 of order unity and screening 
lengths of order of few megaparsecs or smaller. There- 
fore, we can expect that our model predictions differ from 
the ACDM cosmology only on scales ^ lh~~^ Mpc, while 
on larger scales these two models are indistinguishable, 
unlike other modifications of gravity, considered recently, 
for example the DGP model [3l|, f{R) theories , mod- 
ifications of the Newtonian gravity on megaparsec scales 
33, m, [m, or MONDian cosmological simulations 

Here we study only the distribution of the dynami- 
cally dominant dark matter particles. We will study the 
baryon distribution as well in our future work. 



III. SIMULATIONS 



In this section we describe our numerical experiments. 



A. Initial Conditions 

To set up the initial conditions, we have to define the 
power spectrum of the dark matter density fluctuations. 



P{k) - (I'JkP) 



(17) 



We use a power spectrum, derived from the cmbf ast code 
by Seljak & Zaldarriaga [11] with cosmological parame- 
ters h = 0.7, = 0.3, f^A = 0.7 and erg = 0.8. The 
last in this set of parameters is the present value of the 
root-mean-square density contrast of dark matter spatial 
fluctuations within a 8 h"^ Mpc sphere. This is the con- 
ventional normalization parameter and a measure of the 
degree of inhomogeneity of the dark matter distribution. 

The resulting power spectrum, together with the 
PMcode by Klypin & Holtzman [39| is used to displace 
particles from their regular lattice positions, following the 
Zcl'dovich approximation (see Ref. [13] for details). The 
number of individual simulations for each set of model 
parameters has to be large enough to allow a decent aver- 
age over simulation-to-simulation phase fluctuations. We 
decided that for our purposes "a large enough number" is 
10 realizations for 200/i~^ Mpc box simulations. In this 
manner we have obtained an ensemble average of simu- 
lations in the large box. 



B. The codes 

We use the freely available codes: (1) AMIGA {Adap- 
tive Mesh Investigations of Galaxy Assembly) by Knebe, 
Green, Gill & Saar which is the successor of the MLAPM 
code 
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and (2) GADGET2 by Volker Springel [4^, l43|. 
AMIGA is a Particle Mesh code with implementation of 
the Adaptive Mesh Refinements { AMR) technique to ob- 
tain high force resolution. GDAGET2 is a Tree-Particle 
Mesh code. We use AMIGA'S pure Particle Mesh (PM) 



kernel for large box simulation, reducing the simulation 
time at the expense of the force resolution. For study of 
the halo clustering properties we used GDAGET2. To ac- 
commodate for the poor mass resolution we have divided 
our simulations into two sets. 

The first set of simulations was used to study the power 
spectrum and the spatial correlation function of the dark 
matter density fluctuations and density probability dis- 
tribution functions. The simulation box in this series 
of simulations has a width 200/i~^ Mpc, allowing proper 
treatment of the fundamental mode of density pertur- 
bations, which remains in the linear regime at redshift 
z = 0. These are pure PM simulations. 

In the second set of simulations we have used a box of 
width 25/i~^ Mpc. These simulations are used to study 
the cumulative halo mass function and the redshift evo- 
lution of halo abundances and p{5). This approach pro- 
vides a better force and mass resolution, but due to the 
smallness of the box we lack some power in large scales. 
For a discussion of the influence of the box size on the 
dynamics and statistics of simulations, see Ref. ji^-fisl]. 

In each experiment we save the particle positions and 
velocities at redshifts 5; 3; 2; 1; 0.9; 0.8; 0.7; 0.6; 0.5; 
0.4; 0.3; 0.2; 0.15; 0.1; 0.05, and the redshift of the final 
output, z = 0. This archive is used to study the evolution 
of the dark matter distribution with redshift. For an 
experiment, involving Np dark matter particles in a box 
of size L at present, the particle mass, mp, is given by 



^p = (p) {L^Np). 



(18) 



For our simulations, N,, = 128^ for 200 h-^Mpc box 
and Np = 256^ for 25 h'^ Mpc. The other important 
simulation parameters are the force resolution e, the in- 
terparticle separation. 



and the Nyquist wavenumber. 



«Nyq 



tt/£ 



(19) 



(20) 



We list the above parameters, evaluated for the large 
and the small box, in Table [D 



C. The Green's Function 

For conventional gravity in Fourier space, the discrete 
Poisson equation can be written as 



— (>k t/k , 

2a 



(21) 



where is the Green's function. We use the seven-point 
finite-difference approximation to the Green's function to 
solve the Poisson's equation in a cubic box with periodic 
boundary conditions. It is defined for a discrete set of 
arguments. 



q = {q,} = k.L/(27r) 



(22) 
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TABLE L Simulation parameters. L is the box size [/i~^Mpc]; z is the initial redshift; nip is the mass of a single particle 
[h~^ Mq\; e is the force resolution [/i~^kpc]; I is the mean interparticle separation Mpc]; /sNyq is the Nyquist wave 
number [/iMpc~^]; Ts is the screening length [h~^ Mpc]; and (3 is the relative strength of the scalar force. 



L z 


mp 


£ 


I 






P 


200 30 


3.18 • W 


800 


1.563 


2.01 


1;2;5 


-0.5; 0; 0.2; 0.5; 0.7; 1 


25 60 


7.7- 10' 


10 


0.097 


32.1 


0.5; 1; 1.5; 2 


0; 0.2; 0.5; 1 



where L is the present size of the comoving simulation 
box. The subscripts j = 1,2, or 3 denote the three di- 
mensions in k space. The number of grid cells in each 
dimension is M (256 in our simulations), so the cell num- 
bers assume integer values in the range 



1] 



1,2, 



(23) 



Each integer triple defines the grid point, at which the 
Green's function is evaluated: 



AA2^sin2(^q,/AA) 



(24) 



Using the equation. ([9]), we modify Green's function to 
get the proper potential for our model of scalar interac- 
tions: 



-jscalar 



1 



1 + {krsY 



(25) 



The AMIGA code and PM part of the GADGET2 code uses 
equation ()21|) and the fast Fourier transform technique 
to evaluate the gravitational potential at grid points. 
GADGET2 use Oct- Tree algorithm to calculate short dis- 
tance part of the particle forces, we have modified this 
part of the code using equations (|5l6p accordingly. Then 
particle positions and momenta are updated in the stan- 
dard way for N-body algorithms. 



D. Particles and Halos 

To identify collapsed objects from now on called 'ha- 
los' we have used AMIGA'S Halo Finder (AHF) by Knoll- 
man & Knebe Q| which is the successor of the MHF 
halo finder by Gill, Knebe & Gibson The AHF uses 
AMR to find halo centers, then it probes the halo den- 
sity profile around each center in nested radial bins until 
the spatially averaged density contrast reaches the virial 
overdensity, A. At z = 0, in a ACDM universe [48| . 
A = 340. Given our mass resolution in the bigger box, we 
can expect to follow the assembly of halos, correspond- 
ing to clusters and superclusters, and study the statistics 
of clustering at large scales. Simulations in the small 
box should be good enough to investigate the assembly 
of much less massive objects, like halos of clusters and 
galaxies. 

It is important to bear in mind that at small scales 
we are limited by force resolution. The smallest size of 



gravitationally bound objects that can form in our box 
is 2e (see Table |T] and Ref. [i^). We also suffer from the 
well-known problem of overmerpinq described in great 
detail by Klypin et a/, in Ref . [50j. As a consequence, our 
simulations underestimate low-mass object abundances. 



IV. RESULTS 

In this section we present the results obtained in our 
numerical experiments. We provide maps of the spatial 
distribution of dark matter particles as well as clumps of 
particles, called halos. We also study different statistical 
measures of clustering, such as the power spectrum, the 
two-point correlation function, density probability distri- 
bution function, and the cumulative halo mass function 
with and without the scalar interaction. 



A. The Clustering Pattern 

Figure [T] shows the final particle distribution in 
10h~^ Mpc slices, cut through the centers of 200/i~^ Mpc 
simulation boxes. For clarity each slice shows only 1/10 
of the total number of particles, selected at random [6J| 

The frame, labeled 'ACDM' in Figure [T] assumes (3 = 0, 
while the remaining frames show particle distributions 
for a set of different P parameters and a fixed screening 
length, = 2h^^ Mpc. As a test, we also consider 'anti- 
gravity' with (3 = —0.5. All of the frames in Figure [1] 
have evolved from the same initial state, with identical 
amplitudes and phases of density fluctuations, set up at 
z = 30. 

At first sight, modified gravity simulations look similar 
to the ACDM case. The filaments, voids and high-density 
peaks occupy similar positions in the frames. However, 
close inspection shows that with increasing (3, the voids 
appear increasingly more empty. This phenomenon is 
seen more clearly in Figure^ where we show 1.5/i~^ Mpc 
slices, cut through the centers of our smaller simulation 
boxes (L = 25/i~iMpc). 

Figure [2] also shows that the scalar forces enhance the 
accretion of matter, producing more massive halos at 
small scales. Pancake-likc structures, seen in the (3 = 
frame, are more homogeneous than their counterparts in 
frames with /3 > 0, which show substructure. The frag- 
mentation into subsystems of smaller halos is increasingly 
more pronounced, with increasing values of /3. Quantita- 
tively, we can expect an increase of the amplitude of the 
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FIG. 1: Slices cut through centers of simulation boxes of width 200h~^ Mpc at z — 0. The plots show particle positions, x, in 
the {xi,X2) plane. Each slice has dimensions 200 x 200(/i^^ Mpc)'^ in this plane and a thickness of 10h~^ Mpc along the X3 
coordinate axis. 



power spectrum at small scales, corresponding to wave 
numbers k > r~^, and an opposite effect for f3 = —0.5, 
when the pancakes appear even more homogeneous than 
those generated by Newtonian gravity. 



B. The Power Spectrum 

The power spectrum is a convenient measure of the 
strength of dark matter clustering. It is well constrained 
by redshift surveys of galaxies, such as the SDSS catalog. 
In the longwave tail, corresponding to wave numbers 



O.OlMpc^^/i <k < 0.3Mpc~^/i 



(26) 



this survey provides a reliable estimate of P{k) [51|- In 
this range, the ACDM power spectrum agrees well with 
observations. We also use non-inear fit to the power spec- 
trum presented by Smith et al. in ^52| as a measure of 
theoretical ACDM P(k). This can be used to constrain 
scalar field models. 

In Figure [3] we plot the final power spectra, obtained 
from the simulations with L = 200/i~^ Mpc. The top pair 



of plots shows P{k) and the ratio P/Pm for scalar forces 
with fixed (3 = 0.5 and a varying screening length. The 
pair of plots below was obtained from simulations with a 
fixed screening length, = Ih^^ Mpc, and a varying (3. 

These plots show that the rate of growth of density 
perturbations is more sensitive to the range of the scalar 
field, Tg, than to its strength, as long as /3 > 0. The scalar 
force increases the rate of growth of density fiuctuations 
on comoving scales, smaller than and wave numbers 



k > 



As expected, there is an opposite effect for 



P = —0.5: the rate of growth of density perturbations at 
small scales is reduced. 

We describe the deviations from the Newtonian power 
spectrum at a given wave number by the parameter 



AP 



(27) 



These deviations remain negligible on large scales for all 
of the considered models with one exception: the model 
with Ts — 5ft.~^Mpc. For this model, AP/P is signifi- 
cant on all scales present in the simulation box, including 
the SDSS wave number range. Therefore all scalar field 
models with Ts > 5/i~^Mpc do not appear promising. 
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FIG. 2: Similar results for the smaller box, L = 25h 

To assign statistical significance to their failure to repro- 
duce the real Universe, it is necessary to create a mock 
SDSS survey and reproduce the power spectrum estima- 
tor, used by the observers. We are planning to run appro- 
priate Monte Carlo simulations and address this problem 
in greater detail in a forthcoming paper. 

The model with = 5/i~^Mpc may not be suc- 
cessful in reproducing the real Universe, but it is ex- 
tremely useful in understanding the physics of the scalar 
field because this value of the screening length is much 
larger than the interparticle separation in the big box, 
£ ss 1.6/i~^ Mpc. Therefore, we can resolve the differ- 
ence between purely Newtonian gravity and the scalar 
field. Indeed, for k = l/rg = 0.2 , and /3 = 0.5 , we get 
AP/P = 0.5 (see Figure [3]), and an increasing AP/P 
for larger wave numbers. The gravitational attraction, 
enhanced by the scalar field generates the extra power 
in the plot. To look for a similar jump in amplitude for 
smaller values of , we need a better resolution. So, we 
have to turn to the simulations in the smaller box, where 
£ = 125ft.~^kpc. And we do find a similar feature: for 
P = 0.5 and ^ Ih-^Mpc, we get AP/P = 0.5 at 
fc = = 1/iMpc"^ (see Figure H])! 

The decline of Pn /P with decreasing wave number in 



P = 0.5r..= l.Oh'^Mpc 




P= 1.0r,= 2.0h"^Mpc 




5 10 15 20 25 



Xj [Mpc/h] 

^ Mpc. The thickness of each slice is 1.5h~^ Mpc. 

all considered models is a natural consequence of the 
presence of the Yukawa cutoff, reflected in the equation 

In the opposite limit, with growing wave numbers, all 
power spectra in Figure [3] seem to misbehave. Equa- 
tion suggests that all P^/P curves should flatten for 
large wave numbers, k 3> Instead, we see a decline. 

Since this behavior is in disagreement with gravitational 
dynamics, and since it appears in the range k > /cNyq, 
we can expect that the decline is an artifact of the dis- 
crete nature of the simulation. If this is indeed the case, 
then in the simulation in the smaller box, this artifact 
should move to higher wave numbers. This simulation 
has a smaller interparticle separation and a better force 
resolution. Hence, for wave numbers in the range 

< k < /CNyq , (28) 

we should see a flat section of the P/Pn curve, followed 
by a decline for k > fcNyq . This is exactly what happens 
in Figured) The P^/P ratio rises with k until it reaches 
its maximum, followed by a plateau, and then decline for 
wave numbers above the Nyquist limit. 
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FIG. 3: Power spectra from simulations with L — 
200h~^ Mpc. Vertical dashed lines show the Nyquist wave 
number, while the thin black line is the nonlinear fit for 
ACDM P{k) from Smith et al. We also plot the ratio of 
the scalar interaction induced P{k) to its Newtonian counter- 
part, Pt^{k). For all of the upper pair of plots /3 = 1, while 
the Ts parameter is allowed to vary. For the bottom pair 
rs — lh~^ Mpc is fixed, while the value of (3 changes from 
—0.5 to 1. Each curve was obtained by averaging over 10 
realizations with different initial phases. 



C. The Correlation Function 

Another convenient measure of the strength of cluster- 
ing is the spatial two-point correlation function, 

e(|x-x'U) = (<5(x,i)J(x',i)) . (29) 

It is related to the power spectrum by the Fourier trans- 
form, 

^{x) = (2^)^^ / -P(fc)e*-^d^k. (30) 

Under ideal conditions, when the power spectrum is de- 
termined in the entire wave number range from zero to 
infinity, P and ^ are equivalent to each other. However, 



FIG. 4: Power spectra plots similar to those in figure |3] 
obtained from a set of simulations, using the smaller box, 
L = 25/1-1 Mpc. 

measurements from simulations or galaxy redshift sur- 
veys provide only noisy estimates of P{k) or ^(x) for 
limited ranges of k and x. Therefore, in numerical exper- 
iments, it is safer to estimate ^(x) directly from the par- 
ticle positions in the simulation output, using the expres- 
sion , or its discrete version, based on the probability 
density of finding a pair of particles in a separation range 
from X to X + dx (see Ref. [s^). The results presented 
in this section are therefore not a mere Fourier transform 
of the results discussed in the two previous subsections. 
Following the approach we used to study power spectrum 
deviations from the Newtonian case, here we introduce 
a similar measure of the deviations of the scalar-induced 
correlation function from its Newtonian cousin, ^n- This 
is the quantity 

^ . i^. (31) 

The four frames in Figure [5] show two-point correlation 
functions for DM particles in the 200/i^^Mpc box. The 
top two frames show simulation outputs at z = 0, while 
the bottom frames refer to an earlier epoch, z — 1. For 
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FIG. 5: The correlation function for pairs of DM particles at redshift z — 1 (bottom frames) and z = (top frames). The width 
of the simulation box is 200h~^ Mpc. All plots were averaged over 10 simulations. To avoid overcrowding, we show 1 standard 
deviation error bars for the Newtonian ^{x) only. The error bars for the scalar-induced correlation functions are similar. Their 
size is determined mostly by the number of realizations and by the size of the box. Both these quantities are identical for all 
simulations presented here. 



clarity, we show one standard deviation errors only for 
the Newtonian case. The error bars for the remaining 
models are similar. 

The large amplitude in A<^/^ at separations x < 
Ih^^ Mpc is probably an artifact because the separations 
involved are smaller than i — 1.56h^^ Mpc. Nusser et al. 
(2(| have discovered a similar shoulder in ^(x) in their 
simulations and they provided an interpretation similar 
to ours. They have also pointed out that the correlation 
function in their simulations does not possess a feature 
at X = Vs "despite the scalar force attraction at smaller 
scales." We believe that the source of this problem is 
poor resolution, not dynamics. The range of screening 
lengths they consider is uncomfortably close to their in- 
terparticle separation. Our correlation functions, plotted 
in Figure O suffer from the same problem. As we have 
shown in our discussion of the properties of simulated 
power spectra, this problem can be avoided by improv- 



ing the resolution. For an appropriate choice of the box 
size and the screening length, when > l/ZcNyq, the 
feature at fc = l/rg can be resolved. We therefore have 
no doubt - the feature is there. To rediscover for ^{x) 
what we have already seen for P{k) , we only need higher 
resolution simulations, like those in Ref. .53.] , but with a 
modified Green's function. 

Despite the modest resolution of the present results, we 
believe that qualitatively, the stratification of the corre- 
lation functions with respect to and (3, seen in Figure 
[S] reflect the true dynamics. Note that for (3 < , when 
gravity is weaker than Newtonian, we get A^/^ < 0, 
which is dynamically reasonable. 

Our results also differ from those of Nusser et al. at 
larger separations, x > 3h~^ Mpc, where according to 
their Figure 7, A^/^ < 0. This does not happen in 
our simulations unless /3 < 0. This discrepancy, however, 
becomes statistically insignificant if we use our error bars 
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FIG. 6: The density probability distribution functions obtained from 200h ^ Mpc boxes. Results for fixed Ts = Ih ^ Mpc. Left 
panel shows p{S) calculated using the grid cells width of 5h~^ Mpc, while in the right panel are functions calclated with cells 
width of 2.5/1" ^ Mpc. 



for guidance (Nusser et al. do not provide error bars for 
their plots). 

As we have aheady mentioned, for separations x > r^, 
A^/^ is consistent with zero for all models considered 
here. This is good news, since the observed ^(x) at large 
separations agrees well with the ACDM model. 

Last but not least, it is worth noticing that in Figure 
\5\ Af/^ is higher at z = 1 than at z = 0. The growth 
of Newtonian ^ is retarded with respect to the scalar- 
induced ^ at high redshift. Later, 'catches up' with the 
scalar ^. As a consequence, A^/^ is reduced. A possible 
explanation for this phenomenon is that in the presence 
of scalar forces, rapid matter accretion and formation 
of virialized objects occurs at earlier epochs than in the 
standard gravity model. This interpretation is supported 
by our direct analysis of the process of halo formation, 
presented below. 

D. Emptiness of the voids 

As a measure of the emptiness of the voids in our sim- 
ulations we will use the density contrast probability dis- 
tribution function p{d). To calculate these functions we 
extrapolate particle positions to the uniform lattice grid 
to obtain density. After conversion to the density con- 
trast we compute p{d) numerically in the usual way, by 
using 

dp{d6) — , (32) 

where the sum runs for all the cells, and Nc is the total 
number of cells and Si is a density contrast in a i-th 
cell. We have calculated the density contrast probability 
distribution function for the large and small boxes. The 
results are presented in figures [5] and [T] 



In Figure[S]we show p{5) calculated in the 200/i^^ Mpc 
box. Functions in the left frame was calculated using the 
5h~^ Mpc grid cell width, while on the right frame we 
show density probability distribution functions obtained 
with the 2.5ft.~^Mpc grid cell width. In both cases we 
show the effect of varying the (3 parameter keeping the 
screening length = lh~^ Mpc. In Figure [7] we plot 
probability distribution functions (PDFs) oobtained from 
simulations in the small box. The density field was calcu- 
lated on a uniform grid of cells with length of Ih^^Mpc. 
On the left frame we show functions in models with fixed 
Ts = \h~^ Mpc and varying /3. In the right frame we plot 
the results for models with fixed (3 = 1 and various 
values. 

Clearly we see, that introduction of LRSI stretches the 
p{5) toward a smaller density contrast, while keeping the 
relative value at high density tail close to the ACDM case. 
This is striking confirmation of what was expected. LRSI 
produce emptier voids. The smaller the considered cell 
width the more pronounced the effect. An interesting 
effect seen in our distribution functions is the reduced 
skewness. This is different from standard gravity which 
increases the skewness [s^l- We will study this effect in 
the future. 



E. Halos formation in the Small Box 

In this section we study the impact of modified gravity 
on the halo formation process. To identify halos, we ap- 
ply the halo finder program (the AHF, introduced earlier) 
to output snapshots from the 25h~^ Mpc simulations. We 
study all halos with more than 20 particles, hence our 
minimal halo mass is 1.5 • 10^ h~^MQ. 
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FIG. 7: This figure is analogous tho the previous one. Here we present the p{5) obtained from 25h~^ Mpc boxes. The cell 
width used for calculation was lh~^ Mpc. LRSI models in the left panel have fixed = lh~^ Mpc, on the right panles Ts vary 
and we keep (3 = 1 fixed. 




5 10 15 20 25 5 10 15 20 25 
Xj [Mpc/h] Xj [Mpc/h] 



FIG. 8: Positions of halo centers obtained from a particular simulation output. The size of the box is L = 25h~^ Mpc. To avoid 
overcrowding, we show only the halos with centers in the range < a;3 < 5h~^ Mpc. Halos with masses greater than 5 • W^^Mq 
appear as circles. Their diameters are proportional to their masses. Less massive halos are represented by dots. The top left 
frame shows the standard gravity case. The halo distribution in all of the remaining frames is plotted for = lh~^ Mpc and 
three different values of /3. The circle above the upper frame represents a halo with a mass of 2 • 1O"M0. 
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FIG. 9: Maps of the halo distribution, derived from the same simulations as those shown in Figure[8l but for halos with masses, 
smaller than M = 2 . 10^'^ Mq. Even smaller halos with masses axe represented by dots. Those with 

masses in the range m < M < M appear as circles. 



1. High Mass Halos 

In Figure [S] we plot the present comoving positions 
X = {xi,X2,X3) of halo centers in the {xi,X2) coordi- 
nate plane of the simulation box. To avoid overcrowd- 
ing, we consider only halo centers that satisfy the condi- 
tion < 0:3 < 5/i^^Mpc. Halos with masses exceeding 
M = 5 • lOi2/i~U'/0 are plotted as circles with diameters 
proportional to halo masses. For reference, the size of 
a circle, representing a halo with a mass of 2 . IO^^Mq 
is shown in Figure [51 above the upper left frame. The 
halos with masses Af < M are plotted as dots. It is inter- 
esting to note that as expected, the scalar interactions 
enhance the ability of the larger halos to accrete matter 
and become even larger. This effect is particularly strik- 
ing when we compare the upper left, Newtonian frame, 
with the lower right frame where — lh~^ Mpc and 
(3=1. The halos seen in the rectangle 

5/1"^ Mpc < xi <10/i"^Mpc , (33) 
15/1"^ Mpc < X2 < 25/1-^ Mpc , (34) 



at lower right have already acreted all debris in their 
vicinity at earlier times, while in the Newtonian frame 
the accretion process is still going on. At the same time, 
the small-scale mass redistribution process does not affect 
the large scale clustering: the cosmic web is clearly visible 
in all of the frames. 



2. Low-mass Halos 

In this section we look at the low-mass end of the halo 
population. Our Figures [S] and [5] are complementary 
to each other. Both are derived from the same simu- 
lation, except this time we show only halos with masses 
M < m = 2 • lO^O/i^^Mo . The circles show halos with 
masses 8 ■ lO^/i^^M© < M < m. The remaining ha- 
los, with masses M < 8 ■ 10^h~^MQ, are plotted as 
dots. As before, the diameters of the circles are propor- 
tional to halo masses. The influence of the scalar field 
is particularly well pronounced when we compare frames 
with scalar field switched off (/3 = 0, upper left) and on 
(/3 = 1, lower right). The low abundance of light halos, 
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FIG. 10: The redshift evolution of the cumulative halo mass functions in simulations with L = 25h Mpc . We show the 
ACDM case and two scalar models with (3 — 1 and two different values of . 



seen here, is consistent with the high abundance of heavy 
halos in Figure [5] Both figures show the rapid accretion 
and massive halo formation at high redshift, induced by 
the scalar force. 



3. The Cumulative Mass Function 

For a more quantitative description of the halo forma- 
tion process, we will now introduce the cumulative mass 
function (cmf), defined as the mean number density of 
halos with masses greater than the argument mass. 



(35) 
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FIG. 11: Redshift evolution of the abundance of low-massive halos. Here, we consider only halos with M < 2 • lO^^/i ^Mq The 
bottom panels show the ratio of the halo abundance to the ACDM case. 



where Nh{> M) is the number of halos with masses 
greater than M, identified within the a comoving volume 
L^. Later we will also consider the mean number density 
of halos below a certain mass threshold, n{< M) . The 
sum of these two densities, multiplied by gives the 
mean number of halos of all masses. 

In Figure fTOl we show the redshift evolution of the cmf 
for several models. We have also plotted a theoretical 
predictions for the ACDM cmf from the Press-Schechter 
formalism [sBl and Reed et al. f56| using publicly avail- 
able fitting formulas kindly provided by Darren Reed. 
The comparison of the ACDM cmf with its scalar coun- 
terparts confirms the results obtained by plotting the po- 
sitions of halos with different masses. All scalar mod- 
els show enhanced abundances of high mass halos and 
reduced abundances of low-mass halos at z = 0. For 
rs — 2/i^^Mpc, the low-mass tail of the cmf drops by 
almost an order of magnitude below the Newtonian cmf. 
We can also notice that ACDM halo abundance seems 
to be underdeveloped compared to scalar-induced cmf 's 
at high redshifts. This clearly implies, in our opinion, 
that structure formation process in LRSl models is much 
more efficient at high redshifts compared to the Newto- 
nian case. The z = frame in Figure [TU] shows another 



interesting property of the models we consider here. For 
halo masses in the galactic range, IO^^Mq to IO^^Mq , 
the scalar-induced halo abundances agree with the New- 
tonian predictions. This is an advantage because the 
galactic number densities predicted by the ACDM model 
agree with observations [57| . 

Another interesting phenomenon is the overproduc- 
tion of high mass halos in the cluster mass range, M ~ 
10^^h~^MQ. The cluster abundance at 2; = is rela- 
tively well known from observations, and it provides a 
strong cosmological test. It has been used in the past to 
exclude the once popular Einstcin-de Sitter CDM model 
dl]. To decide how deadly this may become for scalar 
interaction models, we need a larger box to sample the 
cluster population properly and more particles to keep 
a reasonable force resolution. Clusters are rare objects. 
In a 25h~^ Mpc box, at the cluster mass range, we are 
dominated by small-number statistics. The mean dis- 
tance between a pair of ACDM clusters (as well as real 
clusters in galaxy surveys) is ^ 50h^^ Mpc. Nusser et at, 
who used a 50h~^ Mpc box, found only a small excess of 
the scalar-induced cluster halo abundance over the New- 
tonian case. We plan to study this problem, using higher 
resolution simulations in the near future. 




FIG. 12: Redshift evolution ol the abundance ol high-massive halos. Here we consider only halos with M > 2 ■ 10^^ Mq 
The bottom panels show the ratio ol the halo abundance to the ACDM case. The small box inset in the bottom right Irame is 
a magnification ol the < z < 0.3 region in the same Irame. 



Apart from plotting 7i(> M) as a function of M for 
fixed redshifts, it is also interesting to fix the mass and 
see how the cmf evolves with z. In Figures [TT] and [T^ 
we present the redshift evolution of two measures of halo 
abundances, n(M < 2-\Q^^h~'^ Mq), shown in the former 
figure, and n{M > 2 ■ lO^/i^^M©), shown in the latter 
figure. 

Note that for all scalar models, the low-mass halo 
abundances are lower than in the Newtonian case at high 
redshift. This happens because the enhanced gravity at 
small scales speeds up the halo formation process, hence 
also increase the typical halo mass at high redshifts com- 
pared to the ACDM case. The light halos become extinct 
earlier, because they merge with larger mass halos. This 
process is responsible for evacuating the voids more ef- 
fectively than conventional Newtonian forces. 

On the high mass end, we see two interesting effects: 
First, because of faster accretion, the higher mass halos 
reach abundances, close to their final values at relatively 
high redshifts. Later, their number densities change less 
rapidly with redshift than the ACDM halo number den- 
sity. This is particularly prominent in the right frame 
at the top of Figure [T^l Such a picture is consistent 



with the observational evidence for an uneventful recent 
past of galaxies like our Milky Way and other nearby 
galaxies, allowing merger events only at high redshifts 
[5^[63,[6l|. In contrast, for the Newtonian model, we see 
rapid growth of n{M > 2 • 10^^h~^MQ) with decreasing 
redshift, suggesting that mergers continue to the present 
day 

The second effect, particularly pronounced in the right 
frame at the top of Figure [121 is the convergence of all of 
the n{> M) curves at z = Q. Since our mass threshhold 
for this set of cumulative mass functions is in the galactic 
mass range, this convergence is a promising feature of 
scalar models, because the ACDM abundance of galaxies 
agrees with observations . 



V. SUMMARY 

We have performed N-body simulations of large scale 
structure formation in a ACDM background, with a New- 
tonian potential and a Newtonian force law, modified by 
the scalar interaction. We have studied the spatial dis- 
tribution of dark matter particles and halos. We also in- 
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vestigated statistical measures of clustering, such as the 
two-point correlation function, the power spectrum, den- 
sity probability distribution function and the cumulative 
halo mass function. We find that the scalar interaction 
removes debris from cosmic voids more effectively than 
the standard ACDM model. It also suppresses late accre- 
tion and merger activity; halo formation processes move 
to higher redshifts. These findings agree very well with 
earlier work j^. For the first time in the literature, we 
have also shown the effect of scalar forces on the evolution 
of the power spectrum. We have resolved the boundary 
between pure Newtonian dynamics, and enhanced power, 
generated by the scalar interactions. 

In the near future we plan to run higher resolution sim- 
ulations. We will follow the evolution of the baryon spa- 
tial distribution, as well as the dark matter. We will also 
consider the three-point correlation function and bispec- 
trum - statistics of higher order, than those considered 
here. Finally we expect to constrain the scalar field pa- 



rameter range by direct comparisons of model predictions 
with observations. 
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